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Optimization of the inverse planning becomes critical because it follows the invention of intensity modu- 
lated radiotherapy (IMRT) to shorten the previous ”trial-and-error’” treatment process and increase efficiency. 
In this paper, the inverse planning is used to direct aperture optimization in the ARTS (Accurate/Advanced 


Radiotherapy System). 


The objective function was quadratic, both tolerance and dose-volume constraint 


types are supported. The memory efficient conjugate gradient algorithm is used to cope with its large data. 
Furthermore, to fully exploit the solution space, a shortest path sub-procedure is coupled into the whole algo- 
rithm, thus giving further possibility decreasing the objective function. Two clinical cases are tested, indicating 
that the applicability of this algorithm is promising to clinical usage. 
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I. INTRODUCTION 


Radiotherapy as an effective tool to cure cancer has been 
practiced for decades. Its guiding principle is to give 
sufficient uniform dose to a tumor target while sparing the 
normal organs. In the old days, this was not an easy task, and 
the treatment was basically based on the physician’s experi- 
ence. The turning point occurred at the 1980’s when Brahme 
solved a simplified inverse problem analytically [1], and an 
ideal goal could be achieved with an inverse treatment to the 
radiotherapy problem. After that, different models and meth- 
ods were used to investigate this issue. Another milestone is 
the introduction of “direct aperture optimization” [2], which 
simplified the traditional optimization process yet still guar- 
antee a promising result. This work falls into the catego- 
ry of direct aperture optimization. The choices of objective 
functions are not limited. Quadratic deviation is a natural 
choice and was adopted by many researchers. Other pop- 
ular objectives include linear [3], biological-based [4], and 
mix integer. A natural choice to optimize would be heuristic 
methods. Simulated annealing was adopted first in the op- 
timization of radiotherapy by Webb [5], and other methods 
like genetic algorithm, evolutionary algorithm etc. are abun- 
dant. The natural advantages of heuristic methods are their 
ability to handle complex, non-regular objective functions, as 
they do not need gradient information to perform optimiza- 
tion. The shortcoming is that the optimization takes much 
time. While deterministic methods surpass in this aspect, by 
using gradient information of objective functions and careful 
implementation, they usually have a much shorter optimiza- 
tion time. Theories on deterministic optimization are enor- 
mous from which effective methods can be adopted on the 
problem of radiotherapy optimization. Due to the fact that 
many optimization algorithms only work on small-scale or 
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middle-scale problems, they do not perform well on large- 
scale problems, and the computer would not even store all 
the data needed in the algorithms. The most effective method 
is conjugate gradient (CG) method to this kind of problem. 
Only part of the data needs to be stored in optimization, and 
it converges in a determined number of iterations. But still, 
special modifications shall be made to customize the stan- 
dard conjugate gradient algorithm applicable to this problem. 
A sub-procedure based on shortest path problem is embedded 
into the optimization to improve the outcome. It searches the 
new apertures which have the potential to decrease the objec- 
tive function value further more. 


II. MATERIAL AND METHODS 
A. Problem formulation, notation 


In the Fusion Driven sub-critical System team [6-9], a ra- 
diotherapy system called ARTS [10—22](Advanced/Accurate 
Radiotherapy System) is developed. The dose engine, 
MCEFSPB (Monte Carlo finite size pencil beam) [19], gener- 
ates the dose kernel matrix for inverse planning optimization. 
We consider two types of organs, tumor target (T) and nor- 
mal organ (NO). Sample dose are picked to fully represent 
the whole organ. Beam intensity, which is denoted by X, is 
a vector with thousands of components, and dose kernel is 
denoted by D. Beam direction is {B}. Each beam direction 
is decomposed into several small beam fields, which are call 
apertures, and denoted by {A}. Each aperture corresponds to 
an aperture weight, w, and is described by the multi-leaf col- 
limator’s leaf position, Y vector. A quadratic based penalty 
function was adopted, as follows: 
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f(d) = ôtpr(dr — Pr)(dr — Pr)” + 6xopNo(dno — Pxo)(dno — Pro}, (1) 


d = Dx, (2) 


$O), (3) 


(4) 


w, x20. (5) 


ô(x), which is a Heaviside function as an indicator, shows 
wether the sample point should be given to penalty. Meeting 
the dose prescription at the sample point is not penalized, and 
ô(x) = 0, or else, penalized, and 6(x) = 1. p is the radio- 
sensitivity factor of different organs. P represents physician’s 
prescription dose to the organ, and means upper tolerance, 
low tolerance, or mean dose in different context. 

In Eq. (3), all the apertures gives the beam intensity vector 
involved in Eq. (1). The dimension of dose kernel is typically 
around 5 000x10 000, and poses limit to even the computa- 
tional power of a modern PC. 


B. Constraints 


Optimization problems and radiotherapy inverse planning 
involves constraints. A group of normal organs called OAR 
(Organ at Risk) are supposed to be spared as much as possi- 
ble, and various constraints are imposed, and an upper toler- 
ance dose is limited by 


Doarx < Poar. (6) 


If the geometrical configuration is too complex, the con- 
straints would degrade the treatment outcome in tumor tar- 
get, the physicians would like to weaken the constraints on 
OAR’s. This option is called as dose-volume constraint, in- 
cluding the dose d, and the volume percentage p%. For a typ- 
ical dose-volume constraint, the percent volume of the organ 
receiving a dose over d should be under p%. In the objec- 
tive function, the target volume received a dose d should be at 
least %. Thus, the dose-volume constraint is more flexible 
than the tolerance limit. Small fraction of the OAR to be over- 
dosed may lead to a uniformly high dose distribution. 

The nonlinear dose-volume constraint incorporating into 
optimization will complicate the optimization process. Here, 
a straightforward sorting-based treatment of dose-volume 
was used, the underdose or overdose part will be penalized 
after sorting the dose. 


C. Algorithms 


To alleviate the computational burden, a multiphase strat- 
egy is applied. Each beam direction is initialized with one 
aperture, the aperture weights are optimized by the Fletch- 
Reeves conjugate gradient method. When the first phase 
stops, the second phase of new aperture generating starts. 
New apertures are generated by solving a shortest path prob- 
lem, which will be explained in Section ITE, added into the 
aperture pool, and the whole optimization process continues. 
After generating 5 new apertures or so, the satisfactory opti- 
mization process is obtained. Fig. 1 shows flow chart of the 
optimization process. 
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CG optimization of 
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Fig. 1. Flow chart of optimization. 


D. Aperture weight optimization 


Simultaneous optimization of aperture shape and aper- 
ture weight is possible, Webb [5] used simulated annealing. 
Genetic algorithm and evolutionary algorithm are adoptable, 
the optimization takes longer time. 

Coupling the aperture weight and aperture shape is the 
main reason why this problem is hard. 


BA (7) 
x= YY wd). 


Handling them separately shall be a wise option. Fixing y 
makes the problem a standard quadratic programming prob- 
lem, which is well solved. After aperture weight optimiza- 
tion, a shortest path problem is formulated as a means to fur- 
ther explore the solution space, which is confined by the pre- 
vious limitation on aperture weight. 
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Conjugate gradient method is suitable for solving large s- 
parse problem, and saves both memory and time. Normally, 
conjugate gradient method performs very well for a standard 
positive definite matrix based quadratic programming. 


minx’ Hx + Px. (8) 


This problem is not quite standard. The objective function 
f(d) = 5¢(dy — PT)(dr — Pr)” + 6No(dxo — Pxo)(dno — Pno)? 
is a weighted quadratic form. Although the square term is 
positive semi-definite, the value of Heaviside function ahead 
each term is not alike each time. So, the problem is evolving 
with time. Nevertheless, the direction provided by a conju- 
gate gradient algorithm is still useful for the evolving. The 
non-negativity constraint on x is incorporated in the optimiza- 
tion by setting the negative components to zero. The follow- 
ing conjugate gradient algorithm is used in this work: 


g(d) = f(d) + ĉtoldioi — Prot) (Ato 
d= Dx. 


E. Shortest path problem 


A necessary condition of the minimization of the objective 
function is 


of /dx = 0. (10) 


The results obtained from Section II D do not meet the con- 
dition because compromises are made. To further explore the 
solution space, finding the solution potential missed can be 
achieved by adding the appropriate aperture grids into the 
solution. The incorporation of a sub-procedure is based on 
an intuition that the negative partial gradient of the objective 
function to the beamlet grid may lead to a decrease of the 
objective function. Gradient of the objective function to the 
beamlet, 0f/0x;, is calculated, and the shortest path problem 
is used for each row of each beam field. 


min sum(0f/0x), —x € B. (11) 


The connectivity requirement of aperture is directly con- 
structed to the optimization process. Various strategies can 
be used to define what does shortest path mean. Here, a naïve 
strategy chosen is the smallest connected beamlet array has 
the most potential to decrease the objective function.The so- 
lution of this problem consists of the new aperture which will 
be added to the aperture pool for optimization. This sub- 
procedure is an effective means to further explore the solution 
space. 


HI. RESULTS AND DISCUSSION 


The test cases consist of a prostate cancer case and 
a lung cancer case. Both of them have voxel sizes of 
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(1) Set initial x9, yo and tolerance €; 

(2) Calculate Vf(xo), if |V f(xo)| < € , stop; otherwise, go 
to Step 3; 

(3) Let po = —V f (xo), let k = 0; 

(4) Line search for tk, Xk+1 = Xk + tkPr, if xj < 0, set x; = 0; 

(5) Calculate V f(xk+1), if |V f(xx+1)| < €, stop; otherwise, 
go to Step 6; 

(6) Prsi = -V f(Xk+1) + AkPk, Ak = mf g =k+1, go 
to Step 4; 

(7) Stop. 

Both tolerance constraint and dose-volume constraint are 
supported and handled by incorporating them into the objec- 
tive function using Heaviside function. Dose-volume con- 
straint is realized by penalizing the region overdosed or un- 
derdosed after sorting. 


(9) 


(5mm x 5mm x 5mm). The dose kernel was computed 
based on an Elekta linac 6MeV electron beam profile in 
ARTS. The configuration of the computer used is: Intel Core2 
Quad CPU Q9550 @2.83 GHz, 3.50GB memory. 


A. Prostate cancer 


The first clinical case is a prostate cancer. It consists of 35 
CT scans. Organs delineated are tumor target, rectum, blad- 
der, and femoral head. Fig. 2 shows the distribution. Tumor 
target is in the middle of four normal organs. Bladder and 
rectum are especially close to the tumor target. 


Fig. 2. (Color online) Organ distribution of the prostate cancer. 


Constraints of the organs are given in Table 1. 
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TABLE 1. Constraints of prostate cancer 


Organ Prescription (cGy) 

Tumor 5100 < p < 5500, mean dose 5200 
Rectum At most 30%, larger than 2000 
Bladder At most 30%, larger than 2000 


Femoral head left 
Femoral head right 


At most 10%, larger than 3500 
At most 10%, larger than 3500 


The 7 beam fields are used in the case with an equi-spaced 
angle of 50°, the total maximum apertures are 50. The op- 
timization took 67s. The iso-dose contour and DVH (dose- 
volume histogram) graph are given in Fig. 3. 
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Fig. 3. (Color online) ISO-dose graph (a) and DVH graph (b) of the 
prostate cancer. 


Although the tumor target is very irregular, the 100% curve 
encompasses the region very well. The 90% and 80% curves 
expanded the central region a little, but they avoided the rec- 
tum as much as possible: almost all the rectum region was 
avoided by the high dose field, only a tiny bit was exposed 
to it. Other organs also receive appropriate dose levels. In 
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TABLE 2. Constraints of lung cancer 


Organ Prescription (cGy) 

Tumor target 7000 < p < 7500, mean dose 7200 
Lung left At most 30%, larger than 3000 
Lung right At most 30%, larger than 3000 
Heart p < 4500 

Spinal cord p < 4100 


the DVH graph, the curve of tumor target is steep, indicating 
good dose conformity to the tumor target. The left and right 
femoral head meet the prescription. And rectum and blad- 
der were constrained very well. The optimization outcome is 
satisfactory. 


B. Lung cancer 


The clinical case of lung cancer contains 53 CT scans. 
The organs delineated are tumor target, lung, heart and spinal 
cord. The distribution is shown in Fig. 4. The tumor target is 
in the left lung. 


Fig. 4. Organ distribution of the lung cancer. 


Objective and constraints of the organs are given in Table 2. 

There are 9 beam fields in this case with an equi-spaced 
angle of 40°, the total maximum apertures are 50. The opti- 
mization took 376 s due to its larger size. The iso-dose con- 
tour and DVH graph are shown in Fig. 5. 

The 100% iso-dose curve followed the same track as the tu- 
mor target contour, though its shape is not quite regular. The 
90% and 80% curves expanded a little, but avoid the neighbor 
spinal cord completely, which is exactly what intensity mod- 
ulated radiotherapy should achieve. An interesting fact is that 
the 80% iso-dose curve is elongated at the bottom direction 
because the tumor target sits in the bottom region of the left 
lung, and higher beam energy was given in the reverse direc- 
tion resulting in a higher tail dose deposition. Other organs 
also received appropriate dose level. From the DVH graph, a 
steep curve of tumor target is presented, indicating good con- 
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Fig. 5. (Color online) ISO-dose graph (a) and DVH graph(b) of the lung cancer. 


puter science, the two aspects are essentially in conflict to 
each other, thus need tradeoffs. Generally, optimization prob- 
lems involved in radiotherapy are broad in type, large-scale 


> formity, spinal cord, heart and lung left were confined very 
y= well under the dose constraints. Dose of lung right is not as 
y= low as the other three, but it still meets the dose criteria. 


c9 and often hard to solve. This work addressed the issue as 
O a direct aperture optimization using a deterministic method. 
© IV. CONCLUSION Conjugate gradient algorithm adopted in this work enables us 


Ever since inverse planning concept was introduced into 
the radiotherapy field, optimization seemed to be a neces- 
sary issue which needs to be addressed. Finding an accept- 
able treatment plan in an acceptable time became the main 
purpose of the radiotherapy treatment planning optimization. 
Unfortunately, the same as other computation issues in com- 


[1] Brahme A. Optimization of stationary and moving beam radia- 
tion therapy techniques. Radiother Oncol, 1988, 12: 129-140. 
DOI: 10.1016/0167-8140(88)90167-3 

[2] Shepard D, Earl M, Li X, et al. Direct aperture optimization: 

a turnkey solution for step-and-shoot IMRT. Med Phys, 2002, 

29: 1007-1018. DOT: 10.1118/1.1477415 

Romeijn H E, Ahuja R K, Dempsey J F, et al. A new lin- 

ear programming approach to radiation therapy treatment plan- 

ning problems. Oper Res, 2006, 54: 201-216. DOI: 10.1287/0- 

pre.1050.0261 

[4] Wu Q W, David D, Wu Y, et al. Intensity-modulated radio- 

therapy optimization with gEUD-guided dose—volume objec- 

tives. Phys Med Biol, 2003, 48: 279. DOI: 10.1088/0031- 

9155/48/3/301 

Webb S. Optimisation of conformal radiotherapy dose distribu- 

tion by simulated annealing. Phys Med Biol, 1989, 34: 1349. 

DOI: 10.1088/003 1-9 155/34/10/002 

Wu Y C. Design status and development strategy of 

China liquid lithium—lead blankets and related material 

technology. J Nucl Mater, 2007, 367: 1410-1415. DOI: 

10.1016/j.jnucmat.2007.04.031 

[7] Wu Y C. Design analysis of the China dual-functional lithium 
lead (DFLL) test blanket module in ITER. Fusion Eng Des, 
2007, 82: 1893-1903. DOI: 10.1016/j.fusengdes.2007.08.012 


[3 


= 


[5 


=a 


[6 


=e 


to solve large-scale problem efficiently. The result is promis- 
ing, and can be obtained in a relatively short time. The heuris- 
tic part of the main algorithm may be the hidden reason why 
this algorithm performs better than some results of the previ- 
ous literature. Some other forms of objectives and constraints 
can be explored in future, which could potentially have a bet- 
ter performance. 


[8] Qiu L J, Wu Y C, Xiao B J, et al. A low aspect ratio toka- 
mak transmutation system. Nucl Fusion, 2000, 40: 629. DOI: 
10.1088/0029-55 15/40/3 Y/325 

[9] Qiu L J, Xiao B J, Chen Y P, et al. A low aspect ratio toka- 
mak transmutation reactor. Fusion Eng Des, 1998, 41: 437- 
442. DOI: 10.1016/S0920-3796(98)003 12-3 

[10] Wu Y C, Li G L, Tao S X, et al. Research and developmen- 
t of an accurate/advanced radiation therapy system (ARTS). 
Chin J Med Phys, 2005, 22: 683-690,702. (in Chinese) DOI: 
10.3969/j.issn. 1005-202X.2005.06.001 

[11] LiGL, Song G, Wu Y C, et al. A multi-objective hybrid genet- 
ic based optimization for external beam radiation. Plasma Sci. 
Technol, 2006, 8: 234-236. 

[12] Li G L, Wu Y C, Song G, et al. Effect of objective func- 
tion on multi-objective inverse planning of radiation thera- 
py. Nucl Phys Rev, 2006. 23: 233-236. (in Chinese) DOI: 
10.3969/j.issn. 1007-4627.2006.02.034 

[13] Cao R F, Li G L, Song G, et al. An improved fast and eli- 
tist multi-objective genetic algorithm-ANSGA-II for multi- 
objective. Chin J Radiol Med Prot, 2007, 27: 467-470. (in 
Chinese) DOI: 10.3760/cma.j.issn.0254-5098.2007.05.016 

[14] Li GL, Song G and Wu Y C. Study on hybrid multi-objective 
optimization algorithm for inverse treatment planning of radia- 
tion therapy. Nucl Tech, 2007, 30: 222-226. (in Chinese) DOI: 


010502-5 


WANG Jie et al. 


10.332 1/j.issn:0253-3219.2007.03.016 

[15] Wu Y C, Song G, Cao R F, et al. Development of accu- 
rate/advanced radiotherapy treatment planning and quality as- 
surance system (ARTS). Chin Phys C, 2008. 32: 177-182. 

[16] Cao R F, Wu Y C, Pei X, et al. Multi-objective optimization of 
inverse planning for accurate radiotherapy. Chin Phys C, 2011, 
35: 313-318. DOI: 10.1088/1674-1137/35/3/019 

[17] Pei X, Cao R F, Liu H, et al. Beam orientation optimization 
using ant colony optimization in intensity modulated radiation 
therapy. WASET, 2011, 56: 590 — 593. 

[18] Li G, Zheng H Q, Sun G Y, et al. Photon energy spectrum re- 
construction based on Monte Carlo and measured percentage 
depth dose in accurate radiotherapy. Prog Nucl Sci Technol, 


Nucl. Sci. Tech. 26, 010502 (2015) 


2011, 2: 160-164. 

[19] Zheng H Q, Sun G Y, Li G, et al. Photon dose calculation 
method based on Monte Carlo finite-size pencil beam model in 
accurate radiotherapy. Commun Comp Phys, 2013, 14: 1415- 
1422. DOI: 10.4208/cicp.221212.100413a 

[20] Wu Y C, Xie Z S and Fischer U. A discrete ordinates nodal 
method for one-dimensional neutron transport calculation in 
curvilinear geometries. Nucl Sci Eng, 1999, 133: 350-357. 

[21] Wu Y C. Conceptual design activities of FDS series fusion 
power plants in China. Fusion Eng Des, 2006, 81: 2713-2718. 
DOI: 10.1016/j.fusengdes.2006.07.068 

[22] Wu Y C. CAD-based interface programs for fusion neutron 
transport simulation. Fusion Eng Des, 2009, 84: 1987-1992. 
DOI: 10.1016/j.fusengdes.2008.12.041 


010502-6 


